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It is now recognized that prediction of radiative heating of entering space craft 
requires explicit treatment of the radiation field from the infrared (IR] to the 
vacuum ultra violet [VUV]. While at low temperatures and longer wavelengths, 
molecular radiation is well described by bound-bound transitions, in the short 
wavelength, high temperature regime, bound-free transitions can play an important 
role. In this work we describe first principles calculations we have carried out for 
bound-bound and bound-free transitions in N 2 , O 2 , C 2 , CO, CN, NO, and N 2 + . 

Compared to bound-bound transitions, bound-free transitions have several 
particularities that make them different to deal with. These include more 
complicated line shapes and a dependence of emission intensity on both bound state 
diatomic and atomic concentrations. . These will be discussed in detail below. 

The general procedure we used was the same for all species. The first step is to 
generate potential energy curves, transition moments, and coupling matrix elements 
by carrying out ab initio electronic structure calculations. These calculations are 
expensive, and thus approximations need to be made in order to make the 
calculations tractable. The only practical method we have to carry out these 
calculations is the internally contracted multi-reference configuration interaction 
(icMRCI] method as implemented in the program suite Molpro. 1 This is a widely 
used method for these kinds of calculations, and is capable of generating very 
accurate results. With this method, we must first of choose which electrons to 
correlate, the one-electron basis to use, and then how to generate the molecular 
orbitals. 


In all calculations, we only correlate the valence electrons, i.e. the lsjike orbitals are 
kept doubly occupied in all configurations. We initially considered correlating all 
electrons, which would yield more accurate results, however the extra expense was 
deemed un-warranted, given the difficulties with the icMRCI method for highly 
excited states. This will be discussed in more detail later. 

For all calculations, the one-electron basis set used was based on the Dunning 2 cc- 
pVTZ basis set, which is the smallest one-electron basis set capable of yielding 
reasonable results. We also used the second order Douglas-Kroll-Hess 3 
approximation to treat scalar special relativity. To help describe the outer-most 
electron clouds, the cc-pVTZ basis set was augmented by one extra diffuse function 
for each shell for C and N, and two extra diffuse functions for each shell for 0. The 
extra shells on 0 are required due to the importance of O' configurations. 4 
Furthermore, systems expected to have low-lying Rydberg states (the neutrals, but 
not the cations] had an additional diffuse function added for each shell. 

The molecular orbitals were generated using the following composite procedure. 

We first generated valance molecular orbitals by carrying out a series of complete 
active space (CAS] calculations starting at large internuclear separations (100 a 0 for 
neutrals and 1000 a 0 for ions] and working our way in to small inter-nuclear 
separations. We will use the symbol r for inter-nuclear separation. About 50 



internuclear separations were considered. In these calculations, the 0 2 sjike 
orbitals were kept doubly occupied in all configurations. In the CAS calculations, 
several roots from each symmetry were optimized, with the weight for each root 
computed from the energy difference as described by Deskevich et al. s We based our 
selection of roots on the atomic states that we could reliable converge to. Thus our 
orbitals transform as Coo V for heternuclear diatomics and Dooh for homonuclear 
diatomics, and furthermore lead to the proper atomic degeneracies in the 
asymptotic region. We just considered a single overall electron spin: doublets for 
CN, N2 + , C 0 + , triplets for N2, O2, C2, CO, and quartets for NO. For N2, C2, CO and NO, 
This choice of spin is not the same as the ground electronic state, but this choice 
gives better coverage of possible atomic states, thus it gives better overall results for 
all spins in the icMRCI calculations. The results from one internuclear separation 
were used as initial guess for the next smaller internuclear separation, and the diab 
procedure in Molpro was used to make the molecular orbitals, both active and 
virtual, look as similar as possible as the internuclear separation changed. 






Once the valence molecular orbitals were determined, we determined low-lying 
Rydberg orbitals for the neutrals. For these calculations we started at the smallest 
inter-nuclear separations and worked our way to larger r values . The valence 
molecular orbitals were kept frozen at their previously determined values, and in 
the CAS calculations we restricted the configurations so that exactly one electron 






was kept in the Rydberg orbital of interest. For N2, O2, C2, and CO, we used the 
lowest lying Rydberg singlet states, while for CN and NO we used the lowest lying 
Rydberg doublet states. In a sequential manner, we generated two a Rydberg 
orbitals and both components of the n Rydberg orbital. As with the valence orbitals, 
we used the diab procedure in Molpro to ensure that the Rydberg orbitals smoothly 
changed with inter-nuclear distance. 

Now we turn to evaluating the energies and coupling matrix elements. We carried 




out icMRCI calculations for using an active space that consisted of all valence 
orbitals and all Rydberg orbitals: however 0 2s valence orbitals where kept doubly 
occupied in the reference configurations and a maximum of one electron was 
allowed on the Rydberg orbitals. For these calculations, we used a modified version 
of Molpro that allowed the selection of A for the state of interest. Thus* X states were 
computed from separate icMRCI calculations then A states, etc. Nonetheless, it is 
necessary to extract multiple roots from the icMRCI calculations. These calculations 
can be carried out in two ways: the coupled state method, in which all roots are 
computed together and The projected state method, in which all roots are computed 
sequentially. The first approach is the default procedure in Molpro and has the 
advantages that the computed wave functions are orthogonal and have lower 
energies, and hence are more accurate, than the projected state energies. It has the 
disadvantage that the cost of calculation grows significantly as the number of states 
increases. Much more troublesome , however, is the fact that changes in the 
character of higher states can lead to discontinuous changes in lower state energies. 
The advantage of the projected state calculations is that the cost is linear in the 
number of roots, and state energies are continuous as the internuclear separation 
changes. The disadvantages of the projected state calculations are the electronic 
wave functions are no longer orthogonal and* if one is unlucky, a state flipping can 
occur that causes the calculations to fail. Thus neither method will always give 
reliable results. For the most part, we use the projected state method, with the 
coupled state method used selectively for O 2 and CO+. The unpleasant features that 
occasionally occur in high lying roots in the figures shown is mainly due to this 
problem. 
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Using these electronic wave functions, we computed the energies, the overlaps, the 
dipole moment, the quadrapole moment, the magnetic dipole moment, and spin- 
orbit matrix elements for all possible couplings. Furthermore, to ascertain the 
relative phases of matrix elements and to transform from the adiabatic to a diabatic 
basis, we computed the overlaps of the wave functions from adjacent inter-nuclear 
separations. 

Now let make a slight digression. The energies computed as described above are 
known as adiabatic energies, for the nuclei are "clamped" to fixed positions in the 
calculations of the electronic energies. However, for accurate calculations, it is 



important to consider the effect of unclamping the nuclei on the electronic 
properties. This is most readily shown in the 2 Yl manifold of NO: in the figure the 
symbols are the adiabatic energies, and the lines are Molecular diabatic energies. 

The adiabatic energies show numerous avoided crossings as the Rydberg state 
comes down in energy . At these crossings, the electronic wave function changes 
rapidly from one side to the other, and thus the coupling between the internuclear 
coordinate and the electronic wave functions will be very large. Even though this 
coupling can be included in the calculations in a fairly straightforward manner, 
other complications with using the adiabatic electronic functions remain. Namely 
the matrix elements of other operators, such as the dipole moment, have very 
complex behavior in the vicinity of an avoided crossing, which makes reliable 
interpolat ion difficult . It is thus advantageous to transform to a different electronic 
basis that does not have these problems because the character of the electronic 
wave functions has been forced to vary slowly and smoothly. This type of basis is 
called a diabatic basis. Several different choices for the diabatic basis are possible, 
and the particular basis we use we call the Molecular diabatic basis. This diabatic 
basis is formed recursively requiring that the adiabatic and diabatic basis be the 
same at the inter-nuclear separation where the adiabatic potential energy curve has 
its deepest minimum. In practice, we start with the internuclear separation at the 
lowest minimum, then propagate both forward and backwards in r. The propagation 
involves a series of Jacobhtype rotations, applied to the overlap between the 
electronic functions at neighboring geometries to make it as diagonal as possible. 
The Jacobi rotations are only applied to one side of the overlap matrix, and each 
individual rotation angle is chosen to minimize the off diagonal overlap of each pair 
of states. Sweeps through the matrix are carried out until the off diagonals of the 
matrix stop decreasing. In the diabatic basis, offidiagonal coupling due to electronic 
energy is non-zero, but all curves are smooth and easy to interpolate. The 
disadvantage of this basis is that convergence to the adiabatic atomic energies is 
rather slow in the asymptotic region. However this is not expected to be a critical 
part of predicting accurate molecular spectra. 


Another thing to note about this figure is the diabatization procedure not only turns 
the sharply avoided crossings into real crossings, but also changes the large r z 
dependence of the ground state potential. In the large r region, the ground state 
potential does not exhibit any obvious avoided crossings, so that the diabatic curve 
differs from the adiabatic curve in this region might be surprising. However some 
consideration of the situation makes the situation more reasonable. In the large r 
region, the potential energy curves can be labeled by atomic quantum numbers, e.g. 
the lowest adiabatic potential curve has the asymptote of 3 P 0 + 4 S N. In the vicinity 
of the minimum, near r=2 a 0t the wave function will reflect a mixture of atomic 
states. For example, due to the higher electro-negativity of 0 over N, there will be 
significant O' N + character in the wave function. In the Molecular diabatic states, the 
fractional character of the wave function is more or less frozen at its value at the 
minimum, thus the energy at the asymptote will be larger than the lowest adiabatic 
curve. This is exactly what we see. What happen to the 3 P 0 + 4 S N state in the 
molecular diabatic basis? After all, even in the vicinity of the minimum, the molecule 
is not purely O N + . The 0 + N component should still be present. The NO figure 
appears to show a deficiency in the molecular diabatic basis. The loss of the 3 P 0 + 4 S 
N state would mean that one cannot use this basis for dissociation. 
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It is interesting to compare the computed potential energy curves with experimental 
data. In the figures, we have placed labels of electronic states at the value of T e and 
r e taken from Huber and Herzberg . 6 For the most part, the agreement is extremely 
good, graphically. The table above with the comparisons between the computed T e 
and those given by Huber and Herzber g in wave numbers . While the differences are 
a small fraction of T e , the absolute errors are sizable enough to be observed in say a 
spectrum from the EAST facility. Thus we shifted the curves for which data was 
available by the negative of the amount listed in the calc-obs columns. A more 
accurate procedure would be to use Tnn, which is a more directly measured quantity, 
but the other approximations in the present work make this level of refinement 
superfluous. 

Now it should be noted that this comparison and correction is not exact, for the 
values of T e computed are just the differences in minimum energies in the Molecular 
diabatic curves, while the T e from Huber and Herzberg is either the Yoo parameter in 
the Dunham expansion of the experimental energy levels for simple cases, or a 
parameter in a complicated effective Hamiltonian for more complex cases. Thus we 


are comparing apples to oranges. Nonetheless it is expected that this correction will 
be reliable enough to yield useful results. 


We next carried out coupled ro-vibrational-spin-electronic level calculations. To 
proceed we build up basis functions for all coordinates except the inter-nuclear 
separation: 
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with P = ±l, A > 0 and E>0 when A=0. Here is the spin-electronic function, 

with p labeling particular electronic states, S is the total electronic spin quantum 
number, A is the absolute value of quantum number specifying the projection of the 
total electronic angular momentum on the inter-nuclear bond axis, E is the quantum 
number specifying the projection of the total electronic spin angular momentum on 
the inter-nuclear bond axis, and x specify spin-electronic coordinates in the frame 
of reference having the body fixed z axis along the inter-nuclear bond, / and M are 
the quantum numbers for the total angular momentum and the projection of the 
total angular momentum on the space fixed z axis. The body fixed and space fixed 
axis are related by the rotation by the Euler angles ^60, P is the parity and x=N e /2, 
where N e is the number of electrons. The factor <^ p is -1 for E - states, and +1 
otherwise. The ^are the Wigner rotation matrix elements as defined by Edmonds. 7 

In this basis, the full relativistic Hamiltonian operator is diagonal in JMP. The non- 
relativistic part of the Hamiltonian is also diagonal in 5. If all off diagonal coupling 
matrix elements are neglected, then this would be what is called the Hund's case a 
basis. However, we do not make this approximation, and rather include all coupling 
terms. Thus we can accurately treat all possible Hund's cases as well as all 
intermediate cases. 


The matrix elements of the nuclear rotational angular momentum contribution 
R(R+1] to the kinetic energy is given by 
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In this expression, we only neglect the * y term, for it is currently difficult to 
calculate, and furthermore only varies by a few cm 1 across the potential energy 
curve. 

Besides the customary electric dipole transition moment, we will consider higher 
order moments. This is because the LBH a^X N2 transition is dipole forbidden. The 
next order moments are the electric quadrapole moment and the magnetic dipole 
moment. In general, the transition matrix elements in this basis are given by 
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is a 3-j symbol. Now this factorizes into a term that contains all the M 
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Now our final wave -functions will consist of linear combinations of the with 
fixed JMP, but this will not effect this factorization, but will change from the labels 
pSAL to v. Then the rate of stimulated emission (aka Einstein A coefficient] from 
electric dipole radiation is 
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the rate of stimulated emission from magnetic dipole radiation is 
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and the rate of stimulated emission from electric quadrapole radiation is 
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We see the following selection rules: for electric dipole radiation (el], P = -P’ and 
\J - J'\ < 1, while for both magnetic dipole (ml] and electric quadrapole (e2] 
radiation we have P - P f , and additionally for e2, \J - J\ = 2 is possible. With all 
other things being approximately equal in atomic units, ml is down from dipole 
allowed transitions by a factor of 1/c 2 , or about 10 -4 . Compared to ml, e2 is down 
about by co 2 , which about another 10 4 in the infra-red region, but nearing unity in 


the VUV. 


So far we have just discussed our treatment of the angular-electronic degrees of 
freedom. The radial degree of freedom, r, is treated using the finite difference 
boundary value method including full coupling. We perform calculations of each 
value of J and P (the radial wave functions are independent of the M quantum 
number] and we explicitly integrate the transition moments over the radial 
functions, thus we do note need Honl-London factors or the r-centroid 
approximation. Furthermore, by fully including all coupling terms, we explicitly 
obtain nonzero lambda doubling splittings, spin-forbidden transitions, and all 
perturbations. 

When solving for the radial functions, we have two situations to consider. In the 
first, for bound levels, the boundary conditions are homogenous as r ->0 and as 
r -> oo. In the second, the total energy is above the lowest dissociation limit, thus the 
large r boundary conditions are inhomogenous. This complicates things 
considerably. 


For bound levels, the homogeneous boundary conditions make the solution of the 
radial functions a banded eigen-value problem. Thus the energies are fixed by the 
boundary conditions. 

On the other hand, for un-bound, or free levels, the energy is a free parameter, and 
we must employ scattering theory to determine the wavefunction. T hen large r 
boundary conditions take the form 

~ ^{ ex p(-*vK' - ex p( ; vK„'} 

where 77 labels channels (channels are the various asymptotically diagonal basis 
functions], = ^2/j^E ~ £ v)’ w ^ ere the reduced mass for radial motion, E is the 

kinetic energy of the relative motion of the free particles. s 71 is the threshold energy 
for channel t], and S is the scattering matrix, which is a complex symmetric matrix. 
For notational simplicity, we have suppressed the energy label from the radial 
function /and the scattering matrix. Now we do not know a priori the values of the 
scattering matrix, so in practical calculations we solve for the functions f w subject 
to the boundary conditions /^(r^) = 8^,. We then analyze the large r behavior of 
the numerical results, and form linear combinations to satisfy the desired boundary 
conditions. It should be noted that the above expressions are only valid when is 
real, i.e. E>£ n . For other channels, the radial functions decay exponentially at large r. 

Now for free states, how should we sample E ? The obvious answer is to use a 
density of energies so that the scattering matrix can be reliably interpolated over. 
This requires knowledge of the formal behavior of the scattering matrix. There are 
poles in the scattering matrix at discrete complex energies E = ^ + iT n 12, where n 
is a counting index. While complex energies are not physically obtainable, these 
poles have important consequences on the real energy axis. Specifically, at energies 
"far" from , i.e. greater than about away, the phase of the determinant of S 

decreases monatonically with increasing E, but near the phase abruptly 


increases by 2n like 2tan l ^T n I2{$g n -is)]. Physically, these poles correspond to 


"almost" bound states with decay life time fi/T n . Thus we need to concentrate our E 
sampling near the <gr . But once again, we do not a priori know the 

We obtain estimates of the <gr by enclosing the system in a finite box, which turns 
all free states into bound states. Then we can determine energy eigen values, and 
the <gr will be nearby some of these eigen values. By evaluating the energy 
derivative of the phase of the scattering matrix, we can isolate the <gr , and then do a 
refinement to determine better estimates of ^ and T n . 

It would seem that knowing <gr and and the f VT1 >, we would be in a position to Formatted: Lowere d by 4 pt J 

compute bound-free transitions just like bound-bound transitions, but further 
reflection reveals additional details that must be addressed. First of all, how should 
the f w be "normalized"? This is not well determined, since they are not square 

integrable functions. Secondly, which of the 77' should we chose to compute 
transition matrix elements? Another way to look at this last problem is f w 
represents incoming flux in channel 77' and outgoing flux in all open channels. But 
the physical situation should have no incoming flux, since all the population would 
be generated by photon absorption. 

The proper solution of this problem is to directly couple the radiation field to the 
nuclear problem using Green's functions. This leads to a driving term in the 
differential equations for the radial functions, and a solution that looks like 



with / indexing the particular bound state for the transition. The amplitude M ' l is 
complex and depends on energy. Note that these differential equations are no more 

difficult to solve than the scattering equations. Then one finds that if<gr corresponds Formatted: Lowered by 4 pt 


to a very nearly bound state, then in the limit of the coupling of the true bound state 
to the continuum becomes zero, then if E is near <^T , 


Formatted: Lowered by 4 pt 
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That is, the emissivity looks like an Einstein A modulated by a Lorentzian, and that 
the Einstein A can be computed via 


A = ^L\K I 


However this simple result is not true in general. This is because the are complex 
numbers, and thus various poles can cause constructive and destructive 
interference. See the following figure for an example: 



In this figure, the red is Jl , which clearly shows first derivative like features 

rather than Lorentzian features. The green is the result of convolution with the 
Doppler broadening. One sees that the broadening for these sharp features can 
result in much of the positive portion of the feature being eliminated. In our 
calculations, we find that this sort of situation is the rule rather than the exception. 

Thus it is not sufficient to just evaluate the amplitude M J at the , but rather it is Formatted: Lowered 4 pt 

important to have more energies in order to nail down the energy dependence of 
the amplitude in order to perform the Doppler broadening calculation. These 



calculations are in progress, but have not yet been completed. 


An additional complication is that in addition to emission from an almost bound 
state, there will be collision-induced emission at that energy. The emission from the 
almost bound state will be the product of the almost bound state population times 

Z l i I 2 

|M ?7 | . In contrast, the collision-induced emission will be the product of the 

Z \ i I 2 

■ See below for an 

v 

example. The red lines are for molecular emission while the other colors are for 
collision-induced emission. 

Thus we only consider bound-bound transitions in this document. In the following 
figures, we show the optically thin absorption spectra of the diatomics computed at 



T=10,000K, which is roughly appropriate for Earth entry. 

We see that in all cases except for N 2 , the dipole transitions dominate the spectrum. 
We also see significant filling in of the spectra between peaks due to weak 
transitions, including spin-forbidden transitions. 
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